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Abstract. We report here on a recent lattice study of the QCD transition region at finite 
Q_j ■ temperature and zero chemical potential using domain wall fermions (DWF). We also present 

^ ' a parameterization of the QCD equation of state obtained from lattice QCD that is suitable 

^ , for use in hydrodynamics studies of heavy ion collisions. Finally, we show preliminary results 

from a multi-stage hydrodynamics/hadron cascade model of a heavy ion collision, in an 
attempt to understand how well the experimental data {e.g. particle spectra, elliptic flow, 
►^ ■ and HBT radii) can constrain the inputs {e.g. initial temperature, freezeout temperature, 

^s ' shear viscosity, equation of state) of the theoretical model. 

^' 
(^ 

(^ ' 1. Introduction 

^D . The aim of the various high energy heavy ion collision (HIC) programs at experimental 

facilities such as RHIC, SPS, LHC, and FAIR is to study the properties of nuclear matter 
under the extreme conditions of high energy and high density. In particular, at sufficiently 

k> , high energy density, it is predicted that normal hadronic matter will undergo a transition into 

?H I a wholly new state of matter, the quark-gluon plasma (QGP) [l|, y], where the constituent 

quarks and gluons will no longer be confined within hadronic states and chiral symmetry is 
restored. 

In principle, the properties of the QGP can be calculated directly from the underlying 
theory. Quantum Chromodynamics (QCD). In practice, ab initio calculations via lattice QCD 
require very large-scale computing resources. It is only recently that high-precision lattice 
QCD results [3|-la] have become available that have the potential to quantatively constrain 
models of heavy ion collisions. In Sec. [2]we discuss a lattice calculation [9|] that uses the domain 
wall fermion method[lO| to calculate the crossover temperature of QCD at zero chemical 
potential. In section [3] we present a parameterization of a high-precision lattice calculation of 
the QCD Equation of State (EoS) that is useful for modeling heavy ion collisions. 

In addition to theoretical calculations of the QGP, a robust understanding of the dynamics 
of a heavy ion collision is needed in order to translate experimental results into constraints on 




Figure 1. The disconnected chiral susceptibility, Xi/^"^ ^nd the subtracted chiral condensate, 
Ai^s/T^- as a function of the bare lattice coupling /3 = Q/g^- The peak in the disconnected 
chiral susceptibility corresponds to T w 170 MeV. 



the physical properties of hot QCD matter. Several approaches, such as Boltzman transport, 
hydrodynamics, and hadronic cascade |llMl4l| have been used to model heavy ion collisions. 
However, the relevant physics changes so drastically over the lifetime of a HIC that it is difficult 
to capture all the aspects of a collision with a single model. In Sec. H]we present preliminary 
results on a multi-stage model of a HIC that incorporates Glauber initial conditions with pre- 
equilibrium flow, 2-d viscous hydrodynamics, Cooper-Frye freezeout, and a hadronic cascade. 



2. QCD Transition using domain wall fermions 

The location of the QCD crossover has been a subject of much debate in the past several 
years [J-[6| • Because of computational expediency, many of these high-precision studies have 
been done using the staggered fermion formulation. Although computationally inexpensive, 
staggered fermions have the known flaw that they only preserve a U{1) subgroup of the full 
SU{Nf) chiral symmetry at finite lattice spacing. This manifests itself in large lattice artifacts 
in certain quantities, e.g. non-degeneracy in the pion spectrum. Domain Wall Fermions 
(DWF) preserve the full SU{Nf) chiral symmetry at finite lattice spacing [lO|. Thus, a DWF 
study at finite temperature is useful to test the robustness of the recent staggered studies of 
Tc and helps us to better understand the role of chiral symmetry in the QCD crossover. 

In our DWF calculation [9|], we perform simulations at 11 different temperatures in the 
range 120 MeV < T < 220 MeV at zero chemical potential. We use a temporal lattice size of 
Nt = 8, which is related to the temperature via T = l/{Nta). These temperatures correspond 
to lattice spacings in the range 0.11 fm < a < 0.21 fm. 

Various quantities are used to locate the chiral-symmetry restoring and deconfinement 
transitions. For the chiral transition, we use the subtracted chiral condensate. 
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where the subscripts I, s denote the light and strange quarks, respectively. In the phase where 
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Figure 2. The isospin susceptibility, C2 Figure 3. The electric charge susceptibil- 



versus bare coupling f3 = 6/g^ 



ity, C2 versus bare coupling /3 = G/qq. 



chiral symmetry is restored (T > Tc), we expect A; ,,(r) to vanish, while for T < Tc, it 
is non-zero. In addition to the chiral condensate, we also calculate the disconnected chiral 
susceptibility. 
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This quantity should exhibit a peak at T = Tc- Fig. [T] shows both 
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seen, the disconnected chiral susceptibility has a clear peak near /3 = 6/gQ 



As can be 
2.03, where go 
is the bare coupling on the lattice. This corresponds to T ~ 170 MeV. 

For the deconfinement transition, we calculate the electric charge susceptibility (c^) and 
the isospin susceptibility (C2): 
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where fiQ and fij are the electric charge and isospin chemical potentials and In Z is the QCD 
partition function. As can be seen in Figs. [2]and[3l these deconfinement observables show a 
rapid rise in the same general region as those observables related to chiral symmetry. 

Although the chiral susceptibility has a clear peak near T ~ 170 MeV and the electric 
charge and isospin susceptibilities imply Tc in the same region, our calculation still suffers 
from significant uncertainties. The spatial volume used is quite small, from V = (1.7 fm)^ up 
to y = (3.1 fm)^, so our result may be polluted with uncontrolled finite- volume effects. In 
addition, the light quark masses are not physical, but correspond to m.^ = 300 MeV near the 
susceptibility peak. Furthermore, the simulations are not done on a line of constant quark 
mass, so that m.j^ is effectively much heavier at low temperatures than at high temperatures. 
These effects introduce significant uncertainty into our estimate of Tc, so we are only able to 
quote a wide range of possible values: 150 MeV < T < 190 MeV. Unfortunately, this does 
not do much to discriminate between the results of recent lattice calculations, but we hope to 
reduce these uncertainties in future calculations. 



3. Parameterization of the lattice QCD equation of state 

There have been several attempts to parameterize the lattice results for the QCD equation 
of state (EoS) in terms of a continuous function that can be conveniently input into 
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Figure 4. Lattice QCD results of the interaction measure (e — 3p)/T^ with corresponding fit 
to the parameterization I{T) = {l/Iiow{T) + l/Ihigh{T))~^. 



hydrodynamic models. Our parameterization reproduces the expected asymptotic behavior of 
the QCD equation of state in the limits T — t- and T — )• oo, smoothly joining them together 
in the intermediate regime near the QCD crossover. 

In lattice calculations, the interaction measure (I(T) = e — 3p) is the thermodynamic 
quantity that is naturally calculated, from which all other thermodynamics quantities such 
as the pressure p, energy density e, entropy density s, and speed of sound Cs can be derived. 
Thus, it is sufficient to parameterize the interaction measure. 

In the high temperature limit, the interaction measure is expected to have the form 15l.ll6|: 
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The lattice QCD data agrees well with this expected behavior. At low temperatures, QCD 
can be described as a gas of mesons and baryons. A widely-used approximation is to assume 
that this gas is non-interacting. This is called the Hadron Resonance Gas (HRG) model [17|]. 
The existing lattice results give an interaction measure that tends to be lower than that 
given by the HRG model at low temperatures. Thus, in our parameterization, deviations from 
the HRG model at low temperature are included in order to accurately fit the lattice data: 
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The low and high temperature parameterizations are combined via: 
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and the parameters (ai • • • 05) can be extracted through a fit to the lattice data. Note that 
fixing 03 = 1 is also an option, if one wants to force the parameterization onto the HRG result 
as T ^ 0. 
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Figure 5. Comparison of the energy density, e and three times the pressure, 3p, for our 
lattice parameterization, the EoS from Laine and Schroeder [18], and the Hadron Resonance 
Gas Model with a 1st order transition. 
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Figure 6. Comparison of the speed of sound squared, Cg, for our lattice parameterization, 
the EoS from Laine and Schroeder [18|], and the Hadron Resonance Gas Model with a 1st 
order transition. 



Table [T] tabulates the fit parameters (qi • • • 05) from a fit to the lattice EoS for the p4 
action at iV^ = 8 [8|. Fig. U] shows the result of this parameterization plotted along with the 
actual lattice data. Fig. [5] shows a comparison of this lattice parameterization with two other 
forms for the equation of state often used in hydrodynamic modeling - a form proposed by 
Laine and Schroeder [18,], and one where the HRG EoS at low temperature is joined to an 
ideal gas at high temperature with a first-order phase transition at Tc = 165 MeV. Figure |6] 
shows the same comparison for the speed of sound squared, Cg. 
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Table 1. Fit parameters for I(T). 



4. Hydro+Cascade Model 

We introduce a multi-stage model where a heavy ion collision is modeled from the first moment 
of collision until the final state hadrons are effectively free-streaming to the particle detectors. 
In constructing this model, we wish to systematically test the effect of varying various 
model inputs and to see which set of parameters give the best agreement with experimental 
results. Among the things that we wish to test are: the nature of the initial conditions, 
the inclusion of pre-equilibirum flow, the initial temperature Tj, the shear viscosity rj, the 
freezeout temperature Ty, and the equation of state. 

4.I. Initial Conditions 



For initial conditions, we use an optical Glauber model J19l . l20l |. where the initial energy 
density is proportional to the number of participant nuclei, Npart or the number of binary 
collisions, Ncoii- The magnitude of the energy density is chosen so that the maximum energy 
density emax corresponds to some initial temperature, Tj. The impact parameter, b can be 
directly varied to obtain different centralities. 

A commonly used assumption is that the initial flow velocities vanishes, neglecting the 
fact that flow velocities may develop in the first r ~ 1 fm/c, before hydrodynamic evolution 
is applicable [2l|. Recently, it has been found that neglecting initial flow velocities might 
have significant effects on collective observables, particularly those related to the source size 



22,123]. Therefore, we allow the modification of our initial conditions to account for possible 



pre-thermalization flow. 

4-2. Viscous Hydrodynamics 

Because the relativistic generalizations of the Navier-Stokes equation are acausal, it is not 
entirely clear how to formulate relativistic dissipative hydrodynamics. One attempt is the 
Israel-Stewart formalism [26], where a relaxation time is introduced for every transport 
coefficient. However, there is still ambiguity as to which higher order terms to include. In 



the case of vanishing bulk viscosity, however, it has been shown [27|] that the most general 



form includes five additional terms, rn, k, At , Ao, A.s, where k = in flat space-time. This 



formulation was implemented by Romatshke [28|], and is the one that we utilize. 



Romatschke's code [28!] is two-dimensional, where the bulk viscosity {Q) implicitly vanishes, 
but the shear viscosity {rj) can be non-zero. For the relaxation times, we set rn = 6r]/sT, 
Ai = 0, A2 = —2rj/Tn, A3 = 0, which are valid choices for a weakly-coupled plasma. The start 
time for the hydrodynamics evolution is taken to be tq = Ifm/c, with the initial conditions 
discussed above. 

For the equation of state, we plan to test four different variations: 1) A HRG EoS for 
T < Tc with a first order transition, 2) the "Laine-Schroder" EoS [181 ] used in [28], 3) a 
parameterization of the lattice QCD EoS discussed in Sec. [3l and 4) the parmeterization 
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Figure 7. vr^ invariant yield. Exper- 
imental data from PHENIX results in 
the 10-15% centrality bin ^. Model 
results obtained with impact parame- 
ter b = 4.4fm, Ti = SOOMeV, Tf = 
ISOMeV, without initial flow and using 
the "lattice- inspired" EoS. 
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Figure 8. vr^ elliptic flow (^2)- Exper- 
imental data from PHENIX results in 
the - 10% centrality bin ^. Model 
results obtained with impact parame- 
ter b = 4.4fm, Ti = 300MeV, Tf = 
ISOMeV, without initial flow and using 
the " lattice- inspired" EoS. 



of the lattice QCD EoS with 03 = 1 so that it is smoothly joined to the HRG EoS at low 
temperatures. A comparison of the first three of these can be seen in Figs. [5] and [6l 

4-3. Cooper-Frye Freezeout 

In modeling the freezeout of hadrons from the QGP, we use sudden freezeout and the Cooper- 
Frye prescription. In this method, the hadrons are frozen out on a hypersurface in (x, y, z, t) 
of constant temperature T = Tf, or equivalently constant energy density. This freezeout 
temperature can be varied to study the effects of early or late freezeout. The single particle 
spectrum using the Cooper-Frye method is given by: 
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where d is the degeneracy factor, dT,^ represents the normal to the freezeout hypersurface, 
and f{x^,p^) is the phase space density with non-equilibrium corrections [30]. 

Because the hydrodynamic evolution occurs only in two spatial dimensions, we assume 
boost- invariance along the longitudinal direction to produce the full freezeout spectrum for 
the particles. 



4-- 4- Hadron Cascade 

After freezeout, one must take into account final state interactions and feed-down decays 
in order to extract the final particle spectra. In order to do this, we take the particles 
produced at freezeout and evolve them through a hadronic cascade. The code that we choose 
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Figure 9. vr^ Rout- Experimental 
data from PHENIX results in the — 



30% centrality bin [29|]. Model results 
obtained with impact parameter b = 
4.4fm, Ti = 300MeV, T/ = ISOMeV, 
without initial flow and using the EoS 
from Laine and Schroeder. 
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Figure 10. vr^ Rside- Experimental 
data from PHENIX results in the — 
30% centrality bin [29]. Model results 
obtained with impact parameter b = 
4.4fm, Ti = 300MeV, T/ = 150MeV, 
without initial flow and using the EoS 
from Laine and Schroeder. 



is the Ultrarelativistic Quantum Molecular Dynamics (UrQMD) code J13l. Il4l|. a microscopic 
transport code that explicitly takes into account particle decays and hadron-hadron collisions. 



4-5. Preliminary Results 

We have preliminary results for several model runs corresponding to Au+Au collisions at 
y/s = 200 GeV per nucleon pair. The parameters that we have used are summarized in Tab. 
[3 However, we have not yet explored all possible variations of the model inputs. The impact 
parameter of 6 = 4.4 fm. was chosen to match as closely as possible to the mid-centraility 
bins for PHENIX results. 

To match to experimental data, we concentrate on three classes of observables: 1) the 
invariant yield {l/27rpT){dN/dpTdy) at mid-rapidity, 2) the elliptic flow V2, and 3) the HBT 
radii. Rout, Rsidei and Riong- To determine the degree of agreement with experimental, we 
perform a joint fit of both the model and experimental results to a common smoothing 
function. For V2 and the HBT radii, the smoothing function we use is a set of orthonormal 
polynomials, the Chebyshev polynomials. For the invariant yield, we found that the product 
of Chebyshev polynomials with exp {—rriT/T), where mr = m? + Pr, produced the best fits. 
Fig. [7] shows the fit of the invariant yield for vr"*" with one set of model parameters. Fig. [8] 
shows the fit for V2 and Figs. [9] and [10] show HBT radii. 

For each set of model parameters, the joint fit with the experimental data produces a value 
of x^i which is a loose measure of the difference between the model and experiment. Figs. [11] 
and [12] show the x^ distributions for vr"*" for the invariant yield and V2 as a function of Ti and 
r]/s. 

As this is a work in progress, we plan to perform a systematic exploration of the model 




Figure 11. x^ distribution for 

model/data fits in the vr"*" invariant 
yield as a function of Tj and r]/ s. 



Figure 12. x^ distribution for 

model/data fits in tt^ V2 as a function 
of Tj and rj/ s. 



parameter space, and perforin these joint fits with a larger set of experimental results. We 
intend to understand how sensitive the final results are to the various model inputs, and which 
model inputs give the best match to experimental results for various systems, energies, and 
centralities. 



b (fm) 


Npart 


Ti (MeV) 


r]/s 


Tf (MeV) 


Initial Flow 


EoS 


4.4 


270 


250, 300, 350 


0.08, 0.16, 0.24 


150 


yes, no 


Laine & Schroder 



Table 2. A summation of the model parameters that have been explored so far for Au+Au 
collisons at y/s = 200 GeV per nucleon pair. 
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